#region Using Directive

using System;
using System.IO;
using NUnit.Framework;
using dnA.Math;
using dnA.Exceptions;

#endregion Using Directive

namespace dnA.dnAtests.Math {

	// suite of tests for double symmetric alogorithm algorithm
	[TestFixture]
	public class DoubleSymmetricLevinsonTest 
	{

		#region ConstructorNullParameter Test
		
		// Test constructor with a null parameter
		[Test]
		[ExpectedException(typeof(System.ArgumentNullException))]
		public void ConstructorNullParameter()
		{
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(null);
		}
		
		#endregion ConstructorNullParameter Test

		#region ConstructorZeroLengthVector
		
		// Test constructor with a zero length vector parameter
		[Test]
		[ExpectedException(typeof(System.RankException))]
		public void ConstructorZeroLengthVector()
		{
			DoubleVector dv = new DoubleVector(1, 0.0);
			dv.RemoveAt(0);
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(dv);
		}
		
		#endregion ConstructorZeroLengthVector
		
		#region GetVector Test
		
		[Test]
		public void GetVector()
		{
			DoubleVector dv = new DoubleVector(new double[]{1.0, 2.0, 3.0, 4.0});
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(dv);
			DoubleVector dslv = dsl.GetVector();
			Assert.IsTrue(dslv.Equals(dv));
		}
		
		#endregion VectorProperty Test
		
		#region GetMatrix Test
		
		[Test]
		public void GetMatrix()
		{
			DoubleMatrix dm = new DoubleMatrix(3);
			dm[0, 0] = dm[1, 1] = dm[2, 2] = 1.0;
			dm[0, 1] = dm[1, 0] = dm[1, 2] = dm[2, 1] = 2.0;
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(dm.GetRow(0));
			DoubleMatrix dsldm = dsl.GetMatrix();
			Assert.IsTrue(dsldm.Equals(dm));
		}
		
		#endregion GetMatrix Test

		#region Decomposition Test 1

		// test the UDL factorisation
		[Test]
		public void Decomposition1()
		{
			// Toeplitz matrix from Golumb and Loan example 4.7.3
			DoubleVector T = new DoubleVector(3);
			T[0] = 1.0;
			T[1] = 0.5;
			T[2] = 0.2;
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T);

			// expected inverse
			DoubleMatrix X = new DoubleMatrix(3);
			X[2, 2] = X[0, 0] = 75.0/56.0;
			X[2, 1] = X[1, 2] = X[1, 0] = X[0, 1] = -5.0/7.0;
			X[2, 0] = X[0, 2] = 5.0/56.0;
			X[1, 1] = 12.0/7.0;

			// get the factors
			DoubleMatrix U = dsl.GetUpper();
			DoubleMatrix D = dsl.GetDiagonal();
			DoubleMatrix L = dsl.GetLower();

			// calculate inverse by the product UDL
			DoubleMatrix Inverse = U * D * L;

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = System.Math.Abs((Inverse[i, j] - X[i, j]) / X[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15);
		}
		
		#endregion Decomposition Test 1

		#region Decomposition Test 2
		
		// test the UDL factorisation with a bigger matrix
		[Test]
		public void Decomposition2()
		{
			const int N = 100;
			DoubleVector T = new DoubleVector(N);
			for (int i = 0; i < N; i++)
			{
				T[i] = 1.0 / (double) (i+1);
			}

			// calculate inverse by the product UDL
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T);
			DoubleMatrix Inverse = dsl.GetUpper() * dsl.GetDiagonal() * dsl.GetLower();

			// multiply with the original matrix and subtract identity
			DoubleMatrix ZeroMatrix = dsl.GetMatrix() * Inverse - DoubleMatrix.CreateIdentity(N);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(ZeroMatrix[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error (bigger matrix, expect more error)
			Assert.IsTrue(MaxError < 1.0E-12, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion Decomposition Test 2

		#region Decomposition Test 3
		
		// test the UDL factorisation with KMS matrix from Freund Example 3 (error = 1.0E-2, which is bad)
		[Test]
		public void Decomposition3()
		{
			const int N = 121;
			const double e = 1.0E-6;
			
			DoubleVector T = new DoubleVector(N);
			T[0] = e;
			for (int i = 1; i < N; i++)
			{
				T[i] = 1.0 / System.Math.Pow(2.0, (double) i);
			}
			// calculate inverse by the product UDL
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T);
			DoubleMatrix Inverse = dsl.GetUpper() * dsl.GetDiagonal() * dsl.GetLower();

			// multiply with the original matrix and subtract identity
			DoubleMatrix ZeroMatrix = dsl.GetMatrix() * Inverse - DoubleMatrix.CreateIdentity(N);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(ZeroMatrix[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error (this is a very ill-conditioned matrix)
			Assert.IsTrue(MaxError < 1.5E-5, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion Decomposition Test 3
		
		#region Decomposition Test 4

		// test the UDL factorisation with a singular matrix
		[Test]
		public void Decomposition4()
		{
			DoubleVector T = new DoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = 1.0;
			
			// calculate inverse by the product UDL
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T);
			dsl.Compute();
			
			Assert.IsTrue(dsl.IsSingular);
		}
		
		#endregion Decomposition Test 4
		
		#region Determinant Test
		
		// Test the Determinant
		[Test]
		public void Determinant()
		{
			// calculate determinant from diagonal
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(new double[]{1.0, 2.0, 3.0, 4.0, 5.0});
			double dsldet = dsl.GetDeterminant();
			
			// use DoubleMatrix member for alternate calculation which is known to be good
			DoubleMatrix ST = dsl.GetMatrix();
			double det = ST.GetDeterminant();
			
			// check results match
			Assert.IsTrue(System.Math.Abs( (dsldet - det)/det ) < 1.0E-15);
		}
		
		#endregion Determinant Test

		#region SolveVector 1
		
		// solve the linear system from Golub, example on page 197.
		[Test]
		public void SolveVector1()
		{
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(new double[]{1.0, 0.5, 0.2});
			DoubleVector X = new DoubleVector(new double[]{355.0/56.0, -376.0/56.0, 285.0/56.0});
			DoubleVector XCalc = dsl.Solve( new double[]{4.0, -1.0, 3.0});
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = System.Math.Abs((XCalc[i]-X[i])/X[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion SolveVector 1
		
		#region SolveVector 2
		
		// solve the linear system from Bareiss, page 413
		[Test]
		public void SolveVector2()
		{
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(new double[]{1.0, 2.0, 3.0, 4.0});
			DoubleVector X = new DoubleVector(new double[]{1.0, 2.0, 3.0, 4.0});
			DoubleVector XCalc = dsl.Solve( new double[]{30.0, 22.0, 18.0, 20.0});
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = System.Math.Abs((XCalc[i]-X[i])/X[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.5E-12, "Maximum Error = "+ MaxError.ToString()); 
		}
		
		#endregion SolveVector 2

		#region SolveMatrix 1

		// calculate inverse by solving linear equations with identity RHS
		[Test]
		public void SolveMatrix1()
		{
			// Toeplitz matrix from Golumb and Loan example 4.7.3
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(new double[]{1.0, 0.5, 0.2});

			// expected inverse
			DoubleMatrix X = new DoubleMatrix(3);
			X[2, 2] = X[0, 0] = 75.0/56.0;
			X[2, 1] = X[1, 2] = X[1, 0] = X[0, 1] = -5.0/7.0;
			X[2, 0] = X[0, 2] = 5.0/56.0;
			X[1, 1] = 12.0/7.0;

			// calculate inverse by solving the linear system
			DoubleMatrix Inverse = dsl.Solve(DoubleMatrix.CreateIdentity(3));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = System.Math.Abs((Inverse[i, j] - X[i, j]) / X[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15);
		}
		
		#endregion SolveMatrix 1
		
		#region SolveMatrix 2
		
		// calculate inverse by solving linear equations with identity RHS
		[Test]
		public void SolveMatrix2()
		{
			const int N = 100;						// the size of the matrix
			DoubleMatrix T = new DoubleMatrix(N);

			// generate test matrices
			double order = (double) N;
			for (int i = 0; i < N; i++)
			{
				double row = (double) i;

				for (int j = 0; j < N; j++)
				{
					T[i, j] = order - System.Math.Abs((row - (double) j));
				}
			}
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T.GetRow(0));

			// calculate the inverse by solving the linear equations
			DoubleMatrix Inverse = dsl.Solve(DoubleMatrix.CreateIdentity(N));

			// multiply toeplitz matrix with inverse and subtract the identity matrix
			DoubleMatrix Y = T*Inverse - DoubleMatrix.CreateIdentity(N);

			// find the maximum error (all elements should be close to zero)
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(Y[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-12);
		}
		
		#endregion SolveMatrix 2
		
		#region GetInverse 1

		// calculate inverse using GetInverse member
		[Test]
		public void GetInverse1()
		{
			// Toeplitz matrix from Golumb and Loan example 4.7.3
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(new double[]{1.0, 0.5, 0.2});

			// expected inverse
			DoubleMatrix X = new DoubleMatrix(3);
			X[2, 2] = X[0, 0] = 75.0/56.0;
			X[2, 1] = X[1, 2] = X[1, 0] = X[0, 1] = -5.0/7.0;
			X[2, 0] = X[0, 2] = 5.0/56.0;
			X[1, 1] = 12.0/7.0;

			// calculate the inverse using Trench's algorithm
			DoubleMatrix Inverse = dsl.GetInverse();

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = System.Math.Abs((Inverse[i, j] - X[i, j]) / X[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15);
		}
		
		#endregion GetInverse 1
		
		#region GetInverse 2
		
		// calculate inverse using GetInverse member
		[Test]
		public void GetInverse2()
		{
			const int N = 100;						// the size of the matrix
			DoubleMatrix T = new DoubleMatrix(N);

			// generate test matrices
			double order = (double) N;
			for (int i = 0; i < N; i++)
			{
				double row = (double) i;

				for (int j = 0; j < N; j++)
				{
					T[i, j] = order - System.Math.Abs((row - (double) j));
				}
			}
			DoubleSymmetricLevinson	dsl = new DoubleSymmetricLevinson(T.GetRow(0));

			// calculate the inverse using Trench's algorithm
			DoubleMatrix Inverse = dsl.GetInverse();

			// multiply toeplitz matrix with inverse and subtract the identity matrix
			DoubleMatrix Y = T*Inverse - DoubleMatrix.CreateIdentity(N);

			// find the maximum error (all elements should be close to zero)
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(Y[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-12);
		}
		
		#endregion GetInverse 2

		#region StaticSolveVectorNullParameter 1 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveVectorNullParameter1()
		{
			DoubleVector X = DoubleSymmetricLevinson.Solve(null, new double[]{1.0, 1.0});
		}

		#endregion StaticSolveVectorNullParameter 1 Test
		
		#region StaticSolveVectorNullParameter 2 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveVectorNullParameter2()
		{
			DoubleVector X = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0}, (DoubleVector) null);
		}

		#endregion StaticSolveVectorNullParameter 2 Test
		
		#region StaticSolveVectorDimension Test

		// test mismatching dimensions
		[Test]
		[ExpectedException(typeof(RankException))]
		public void StaticSolveMatrixDimension()
		{
			DoubleVector X = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0, 1.0}, new double[]{1.0, 1.0});
		}
		
		#endregion StaticSolveVectorDimension Test
		
		#region StaticSolveVectorSingular
		
		// test with Toeplitz matrix which has a singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticSolveVectorSingular()
		{
			DoubleVector X = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0, 1.0}, new double[]{1.0, 1.0, 1.0});
		}
		
		#endregion StaticSolveVectorSingular
		
		#region StaticSolveVectorExample 1 Test
		
		// solve the linear system from Golub, example on page 197.
		[Test]
		public void StaticSolveVectorExample1()
		{
			DoubleVector X = new DoubleVector(new double[]{355.0/56.0, -376.0/56.0, 285.0/56.0});
			DoubleVector XCalc = DoubleSymmetricLevinson.Solve(new double[]{1.0, 0.5, 0.2}, new double[]{4.0, -1.0, 3.0});
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = System.Math.Abs((XCalc[i]-X[i])/X[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveVectorExample 1 Test
		
		#region StaticSolveVectorExample 2 Test
		
		// solve the linear system from Bareiss, page 413
		[Test]
		public void StaticSolveVectorExample2()
		{
			DoubleVector X = new DoubleVector(new double[]{1.0, 2.0, 3.0, 4.0});
			DoubleVector XCalc = DoubleSymmetricLevinson.Solve(new double[]{1.0, 2.0, 3.0, 4.0}, new double[]{30.0, 22.0, 18.0, 20.0});
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = System.Math.Abs((XCalc[i]-X[i])/X[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 2.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveVectorExample 2 Test

		#region StaticSolveMatrixNullParameter 1 Test

		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveMatrixNullParameter1()
		{
			DoubleMatrix Y = new DoubleMatrix(2);
			Y[0, 0] = Y[1, 1] = 2.0;
			Y[0, 1] = Y[1, 0] = 1.0;
			DoubleMatrix X = DoubleSymmetricLevinson.Solve(null, Y);
		}
		
		#endregion StaticSolveMatrixNullParameter 1 Test
		
		#region StaticSolveMatrixNullParameter 2 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveMatrixNullParameter2()
		{
			DoubleMatrix Y = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0}, (DoubleMatrix) null);
		}
		
		#endregion StaticSolveMatrixNullParameter 2 Test
		
		#region StaticSolveMatrixDimension 1 Test

		// test mismatching dimensions
		[Test]
		[ExpectedException(typeof(RankException))]
		public void StaticSolveMatrixDimension1()
		{
			DoubleMatrix Y = new DoubleMatrix(2);
			Y[0, 0] = Y[1, 1] = 2.0;
			Y[0, 1] = Y[1, 0] = 1.0;
			
			DoubleMatrix X = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0, 1.0}, Y);
		}
		
		#endregion StaticSolveMatrixDimension 1 Test
		
		#region StaticSolveMatrixSingular

		// test with Toeplitz matrix which has a singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticSolveMatrixSingular()
		{
			DoubleMatrix X = DoubleSymmetricLevinson.Solve(new double[]{1.0, 1.0, 1.0}, DoubleMatrix.CreateIdentity(3));
		}
		
		#endregion StaticSolveMatrixSingular
		
		#region StaticSolveMatrixExample 1 Test

		// calculate inverse by solving linear system with identity RHS
		[Test]
		public void StaticSolveMatrixExample1()
		{
			// expected inverse
			DoubleMatrix X = new DoubleMatrix(3);
			X[2, 2] = X[0, 0] = 75.0/56.0;
			X[2, 1] = X[1, 2] = X[1, 0] = X[0, 1] = -5.0/7.0;
			X[2, 0] = X[0, 2] = 5.0/56.0;
			X[1, 1] = 12.0/7.0;

			// calculate inverse by solving the linear system
			DoubleMatrix Inverse = DoubleSymmetricLevinson.Solve(new double[]{1.0, 0.5, 0.2}, DoubleMatrix.CreateIdentity(3));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = System.Math.Abs((Inverse[i, j] - X[i, j]) / X[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15);
		}
		
		#endregion StaticSolveMatrixExample 1 Test
		
		#region StaticSolveMatrixExample 2 Test
		
		// calculate inverse by solving linear system with identity RHS
		[Test]
		public void StaticSolveMatrixExample2()
		{
			const int N = 100;						// the size of the matrix
			DoubleMatrix T = new DoubleMatrix(N);

			// generate test matrices
			double order = (double) N;
			for (int i = 0; i < N; i++)
			{
				double row = (double) i;

				for (int j = 0; j < N; j++)
				{
					T[i, j] = order - System.Math.Abs((row - (double) j));
				}
			}

			// calculate inverse by solving the linear system
			DoubleMatrix Inverse = DoubleSymmetricLevinson.Solve(T.GetRow(0), DoubleMatrix.CreateIdentity(N));

			// multiply toeplitz matrix with inverse and subtract the identity matrix
			DoubleMatrix Y = T*Inverse - DoubleMatrix.CreateIdentity(N);

			// find the maximum error (all elements should be close to zero)
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(Y[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-12);
		}
		
		#endregion StaticSolveMatrixExample 2 Test

		#region StaticYuleWalkerNullParameter Test
		
		// test Yule-Walker with a null reference
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticYuleWalkerNullParameter()
		{
			DoubleVector Y = DoubleSymmetricLevinson.YuleWalker(null);
		}
		
		#endregion StaticYuleWalkerNullParameter Test
		
		#region StaticYuleWalkerDimension Test

		// test Yule-Walker with order 1 matrix
		[Test]
		[ExpectedException(typeof(ArgumentOutOfRangeException))]
		public void StaticYuleWalkerDimension()
		{
			DoubleVector Y = DoubleSymmetricLevinson.YuleWalker(new double[]{1.0});
		}
		
		#endregion StaticYuleWalkerDimension Test
		
		#region StaticYuleWalkerSingular
		
		// test Yule-Walker with matrix with singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticYuleWalkerSingular()
		{
			DoubleVector Y = DoubleSymmetricLevinson.YuleWalker(new double[]{1.0, 1.0, 1.0});
		}
		
		#endregion StaticYuleWalkerSingular
		
		#region StaticYuleWalkerExample 1 Test
		
		// system values from Golub and Loan, Example 4.7.1
		[Test]
		public void StaticYuleWalkerExample1()
		{
			DoubleVector A = DoubleSymmetricLevinson.YuleWalker(new double[]{1.0, 0.5, 0.2, 0.1});
			DoubleVector X = new DoubleVector(3);
			X[0] = -75.0/140.0;
			X[1] = 12.0/140.0;
			X[2] = -5.0/140.0;

			// find the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = System.Math.Abs((A[i]-X[i])/X[i]);
				if (E > MaxError)
				{
                	MaxError = E;
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15);
		}

		#endregion StaticYuleWalkerExample 1 Test
		
		#region StaticYuleWalkerExample 2 Test
		
		// system values from NAG F04MEF
		[Test]
		public void StaticYuleWalkerExample2()
		{
			DoubleVector A = DoubleSymmetricLevinson.YuleWalker(new double[]{4.0, 3.0, 2.0, 1.0, 0.0});
			DoubleVector X = new DoubleVector(4);
			X[0] = -0.8;
			X[1] = 0.0;
			X[2] = 0.0;
			X[3] = 0.2;

			// find the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = System.Math.Abs(A[i]-X[i]);
				if (E > MaxError)
				{
                	MaxError = E;
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticYuleWalkerExample 2 Test
		
		#region StaticInverseNullParameter Test

		// DoubleSymmetricLevinson.Inverse(DoubleVector T) with null parameter.
		[Test]
		[ExpectedException(typeof(System.ArgumentNullException))]
		public void StaticInverseNullParameter()
		{
			DoubleMatrix Y = DoubleSymmetricLevinson.Inverse(null);
		}

		#endregion StaticInverseNullParameter Test

		#region StaticInverseDimension Test

		// DoubleSymmetricLevinson.Inverse(DoubleVector T) with small dimension
		[Test]
		[ExpectedException(typeof(System.RankException))]
		public void StaticInverseDimension()
		{
			DoubleMatrix Y = DoubleSymmetricLevinson.Inverse(new DoubleVector(new double[]{1.0}));
		}

		#endregion StaticInverseDimension Test

		#region StaticInverseSingular

		// DoubleSymmetricLevinson.Inverse(DoubleVector T) with an ill-conditioned system
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticInverseSingular()
		{

			// setup an ill-conditioned system (second order principal submatrix is singular)
			DoubleVector T = new DoubleVector(3);
			T[0] = 1.0;
			T[1] = 1.0;
			T[2] = 1.0;

			DoubleMatrix Y = DoubleSymmetricLevinson.Inverse(T);

		}

		#endregion StaticInverseSingular

		#region StaticInverseGolumbExample Test

		// DoubleSymmetricLevinson.Inverse(DoubleVector T) with Golumb and Loan example 4.7.3
		[Test]
		public void StaticInverseGolumbExample()
		{
			// Toeplitz matrix from Golumb and Loan example 4.7.3
			DoubleVector T = new DoubleVector(3);
			T[0] = 1.0;
			T[1] = 0.5;
			T[2] = 0.2;

			// expected inverse
			DoubleMatrix X = new DoubleMatrix(3);
			X[2, 2] = X[0, 0] = 75.0/56.0;
			X[2, 1] = X[1, 2] = X[1, 0] = X[0, 1] = -5.0/7.0;
			X[2, 0] = X[0, 2] = 5.0/56.0;
			X[1, 1] = 12.0/7.0;

			// calculate the inverse
			DoubleMatrix Y = DoubleSymmetricLevinson.Inverse(T);

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = System.Math.Abs((Y[i, j] - X[i, j]) / X[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15);

		}

		#endregion StaticInverseGolumbExample Test

		#region StaticInverseBig Test

		// check a large system for accuracy
		[Test]
		public void StaticInverseBig()
		{

			const int N = 100;						// the size of the matrix
			DoubleMatrix T = new DoubleMatrix(N);

			// generate test matrices
			double order = (double) N;
			for (int i = 0; i < N; i++)
			{
				double row = (double) i;

				for (int j = 0; j < N; j++)
				{
					T[i, j] = order - System.Math.Abs((row - (double) j));
				}
			}

			// calculate the inverse
			DoubleMatrix Y = DoubleSymmetricLevinson.Inverse(T.GetRow(0));

			// multiply toeplitz matrix with inverse and subtract the identity matrix
			Y = T*Y - DoubleMatrix.CreateIdentity(N);

			// find the maximum error (all elements should be close to zero)
            double MaxError = 0.0;
            for (int i = 0; i < N; i++)
            {
            	for (int j = 0; i < N; i++)
				{
                	double E = System.Math.Abs(Y[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-12);

		}
		
		#endregion StaticInverseBig Test
		
	}

}

